Directed avalanche processes with underlying interface dynamics 
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We describe a directed avalanche model; a slowly unloading sandbox driven by lowering a retaining 
wall. The directness of the dynamics allows us to interpret the stable sand surfaces as world sheets 
of fluctuating interfaces in one lower dimension. In our specific case, the interface growth dynamics 
belongs to the Kardar-Parisi-Zhang (KPZ) universality class. We formulate relations between the 
critical exponents of the various avalanche distributions and those of the roughness of the growing 
interface. The nonlinear nature of the underlying KPZ dynamics provides a nontrivial test of 
such generic exponent relations. The numerical values of the avalanche exponents are close to the 
conventional KPZ values, but differ sufficiently to warrant a detailed study of whether avalanche 
correlated Monte Carlo sampling changes the scaling exponents of KPZ interfaces. We demonstrate 
that the exponents remain unchanged, but that the traces left on the surface by previous avalanches 
give rise to unusually strong finite-size corrections to scaling. This type of slow convergence seems 
intrinsic to avalanche dynamics. 

PACS numbers: 45.70.Ht, 05.65.+b, 05.70.Np, 47.54.+r 
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I. INTRODUCTION 

Avalanche phenomena are common in nature. Exam- 
ples range from accumulating snow on mountain slopes, 
slow shearing between continental plates [Q, rerouting 
in river networks, to creeping magnetic flux lines in 
super conductors 0|. Following the work by Bak et 
al. |3|, physicists aim to capture the essential aspects 
of such dynamical systems with simple automaton pro- 
cesses, commonly referred to as sandpile models and self- 
organized criticality (SOC). Impressive successes have 
been achieved, like reproducing power-law distributions 
in avalanche events similar to those observed in nature, 
and the start of a classification scheme of such processes 
in terms of so-called universality classes M. Unfortu- 
nately most of these are numerical in nature. Analytical 
exact results remain rare. 

Directed avalanche phenomena form a subclass of these 
SOC processes. Dhar and Ramaswamy introduced the 
first directed sandpile model and solved it exactly ||. 
This was possible because in their model the avalanche 
propagation is governed solely by its two edges, and those 
two follow independent random walk dynamics. Tadic 
and Dhar Q introduced a directed model in which par- 
ticles are allowed to pile up beyond the critical height, 
by replacing the automaton's deterministic toppling rule 
by a stochastic one ||. The density of critical sites 
tunes itself and at distances far from the driving edge 
the propagation of active sites approaches the directed 
percolation (DP) 0] threshold. The scaling properties of 
the avalanche distributions are thus linked to the crit- 
ical exponents characterizing the DP universality class. 
Another example of a stochastic directed avalanche pro- 
cess is the model introduced and studied numerically by 
Pastor-Satorras and Vespignani ||. Similar as in the 
above model by Dhar and Ramaswamy, the stable land- 
scape configurations (between avalanche events) lack in- 
ternal correlations in the stationary state. This allowed 



Paczuski and Bassler 11 and also Kloster et al. jT(J to link 
this dynamic process to so-called Edwards- Wilkinson |ll[] 
(EW) interface growth and to derive the exact scaling ex- 
ponents of the avalanche distributions. 

This novel world sheet type connection between 
avalanche dynamics and interface growth is particularly 
promising, because interface dynamic processes like EW 
and Kardar-Parisi-Zhang [fijj (KPZ) growth are very 
well understood, in particular in 1+1 dimensions (1+1D) 
where the scaling properties are known exactly. However, 
the above models that are linked to EW type growth are 
rather poor examples, because EW growth is described 
by a simple linear stochastic (diffusion type) Langevin 
equation; correlations factorize, and important caveats 
in the relation to avalanche dynamics can be obscured 
by this simplicity. 

We set out to generalize this approach to nonlinear 
interface dynamic processes, and recently introduced a 
directed unloading sandbox model |]l3| in which the two 
dimensional (2D) avalanche dynamics relates to 1+1D 
KPZ type interface growth. We derived exponent rela- 
tions between the avalanche and interface growth scal- 
ing properties, which are generic, and valid beyond our 
specific model. Our numerical results for the avalanche 
distributions (for length, width, depth, and mass) follow 
indeed these exponent relations. Moreover, the avalanche 
critical exponents obey the predicted KPZ values within 
a few percent, an accuracy typical to avalanche simula- 
tions. However, our numerical accuracy is better than 
that; mostly because of a careful finite-size scaling (FSS) 
analysis. The exponents seem to converge to values that 
are slightly different from the KPZ values. 

This left us with a puzzle. What is the origin of these 
small deviations? Is this a fundamental effect; or do the 
exponents ultimately converge to the KPZ values, but 
with unusually large corrections to scaling. In this pa- 
per we address these issues. We also provide a more de- 
tailed discussion of these world-sheet-type relationships 
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between avalanche and interface growth dynamics. Our 
first paper was short and did not include many of the 
details that are crucial for the analysis presented here. 

The fundamental difference between conventional KPZ 
interface growth and avalanche dynamics arises from the 
averaging process over KPZ type space-time world sheets. 
In normal Monte Carlo (MC) simulations of interface 
growth the distribution functions are determined in terms 
of ensemble averages over a set of totally uncorrelated 
space-time MC runs. In contrast, the avalanche dynam- 
ics gives rise to KPZ world sheets that are strongly cor- 
related. Two subsequent MC runs are identical except 
inside a single avalanche area. This difference in averag- 
ing, uncorrelated versus avalanche correlated MC runs, 
therefore emerges as a key issue for understanding the 
scaling properties of avalanche dynamics. This issue did 
not arise in the earlier EW type avalanche models due 
to the linear nature of the EW process. However for 
nonlinear dynamics, like KPZ, avalanche-correlated-type 
sampling could well lead to novel interface scaling expo- 
nents. 

Speaking against a shift in the values of the exponents, 
are arguments like: the KPZ stationary state, i.e., the 
sand surface profile far way from the driving edge, can not 
be affected by the avalanche-correlated-type averaging, 
because large avalanches that span the entire width of 
the box occur periodically. These completely refresh the 
surface far way from the driving edge regularly, and thus 
wipe out all correlations between MC runs. This suggests 
that we are only dealing with much larger than usual 
corrections to scaling. The details are more complex than 
this simple argument, but we will establish that indeed 
the exponent values do not change. 

The paper is organized as follows. In the next section 
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we present the unloading sandbox model. In Sec 
we comment on how directed avalanche dynamics can be 
linked to interface growth in one lower dimension. Next, 



in Sec. IV, we show that in the interface growth interpre- 



tation our specific model belongs to the KPZ universality 
class. In Sec. |y| we derive the generic exponent relations 
between interface growth and directed avalanche dynam- 
ics, and in Sec. [v| we test this numerically for our specific 
model. 

In the second half of this paper we address the small 
deviations in the numerical values of the expone nts fr om 
those of conventional KPZ growth. In Sec. VII we 
present numerical results detailing how the traces left on 
the surface profile by previous avalanches influence both 
the avalanche exponents and the interface growth ones. 
These scars in the rough surface enhance the surface 
roughness. We cast this enhanced interface roughness 
in terms of corrections to scaling, and determine what 
value the critical dimension of the corresponding irrele- 
vant operator O sc (in the sense of renormalization theory) 
should have. Next, we identify the geometric meaning of 
O sc , starting with a study of the one dimensional (ID) 
version of our model where a similar phenomenon takes 




FIG. 1: Sandbox with a slowly lowering retaining wall 



simple random walk, and the avalanche correlated sam- 
pling relates to the scaling properties of merging random 
walkers. O sc represents the distribution of avalanche end- 
points in the ID surface, and can be studied directly from 
the rounding of the surface profile near the driving edge. 
In Sec. [X we return to the full 2D case. The scars of pre- 
vious avalanches form lines on the surface. We identify 
O sc with the angle these lines make with respect to the 



place, Sec. VIII. In ID the interface growth process is a 



direction perpendicular to the driving edge, and confirm 
with an analytic argument that the critical dimension of 
O sc is equal to :r sc = — z with z the KPZ dynamic expo- 
nent. Finally, we summarize our results in Sec. [x| 



II. AN UNLOADING SANDBOX 

Imagine a box filled with granular material, as illus- 
trated in Fig. |l|. One of its four retaining walls is slowly 
lowered, such that the sand spills out from that side, and 
thus slowly unloads the box and establishes a sloped sur- 
face. In the quasistatic limit, the wall moves slow enough 
that the unloading events can be described as distinct 
avalanches. The box can be three dimensional, leading 
to 2D avalanche dynamics on a 2D surface, or can be 2D 
(like in a very narrow box) giving rise to ID avalanches 
on a ID surface. 

Inspired by this we consider a so-called solid-on-solid 
model defined on a 2D lattice. Height variables h(r) are 
defined on a square lattice. We will consider two ver- 
sions of the model. In the continuous height version, 
the heights are real numbers. In the discrete model, the 
heights are integers, h(r) = 0, ±1,±2, The former 
corresponds to a continuous material without internal 
structure, but strong cohesion up to a specific length scale 
s c , while the latter corresponds to layered material where 
the surface height is quantized. 

The 2D lattice is rotated diagonally such that the prop- 
agation direction of the avalanche is along the diagonal 
direction denoted by y. This is the direction in which 
the avalanche will run. Throughout this paper the coor- 
dinate perpendicular to y will be denoted by x. Figure ^ 
illustrates this geometry. 
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FIG. 2: Lattice structure of sandbox model in 2D 



The configurations are subject to the following stabil- 
ity condition. The column of particles on site r = (x, y) 
is supported by the two columns, r/ = (x — l,y — 1) and 
r r = (x + 1, y — 1) directly below it and is stable when its 
height is less than the minimum of the heights at these 
two supporting sites increased by a fixed amount 

h(r) < min [h(ri), h(r r )] + s c . (1) 

s c is a constant. In the version of our model where the 
heights are continuous variables s c represents the only 
length scale in the h direction and can be set equal to 1 
without loss of generality. Throughout this paper we will 
also set s c = 1 in the discrete h model. 

Consider a stable configuration, after t — 1 avalanches. 
The t-th avalanche is triggered at the highest site r = 
(xj, 0), on the y — driving boundary (or, in the discrete 
height model, by randomly choosing one of the high- 
est sites) and reducing its height by a random amount 
< Vt — s c- This likely creates unstable sites in the 
next y = 1 row. Those are updated by replacing their 
height by an amount equal to the lowest of the two sup- 
porting columns in the previous row and then adding an 
uncorrelated random amount < ?y(r) < s c with uniform 
distribution, as 

h(r) ^miii[/i(r,),/i(JV)]+»7(r). (2) 

This updating continues row by row until all the sites 
are stable again. Only after that the next avalanche is 
started. The toppling of a site only effects the stability 
of the two sites immediately above it in the next y-row. 
Therefore we can update the system row-by-row in in- 
creasing order of y. 

Direct experimental realizations of this unloading 
sandbox model are not our immediate concern (the fo- 
cus is on establishing a generic theoretical relationship 
between avalanche dynamics and interface growth), but 
we expect that this model is applicable to actual experi- 
mental unloading sandboxes. One of the most important 
issues in this context is the row-by-row nature of the 
toppling rule. This is a crucial feature for our purposes, 
allowing the identification with KPZ interface growth (in 
the next section). In real unloading sandboxes the sand 
removed from row y rolls down hill and likely disturbs 



the already stabilized lower surface levels. Experimental 
realizations can avoid this from happening, e.g., by choos- 
ing very light grains (compared to the cohesion forces). 
Note that our dynamic rule does not allow the build-up 
of any pockets (deeper than s c ) on the surface that might 
trap such downward rolling grains. 

Conservation laws are crucial to avalanche dynamics. 
Unlike most avalanche processes, our model does not con- 
serve mass while the avalanche propagates. That might 
raise the specter of our model not being (self-organized) 
critical. The connection to KPZ growth (an intrinsic crit- 
ical process) dispels this phantom. Moreover, the global 
slope of the surface is preserved during each avalanche 
run, and conservation of steps in the profile plays the 
role analogous to conservation of mass. 

The analysis of the dynamics involves distribu- 
tion functions of various characteristic features of the 
avalanches, The common examples are: length, width, 
depth, and mass. The avalanche length I will be defined 
throughout this paper as the maximum distance y the 
avalanche travels from the driving edge; the width w as 
the maximum departure of the ^-coordinate (perpendic- 
ular to the propagation direction) from the trigger point 
x-coordinate; the depth S as the maximum height change 
the avalanche creates at any of the affected sites; and the 
mass m as the total amount of material removed by the 
avalanche. 



III. AVALANCHES VERSUS EPITAXIAL 
INTERFACE GROWTH 

The focus of this paper is on how the above avalanche 
dynamics relates to interface growth in one lower dimen- 
sion. Each stable sloped surface configuration of a di- 
rected sandpile can be reinterpreted as a world sheet 
(space-time configuration) of an interface in one lower 
spatial dimension. The direction in which the avalanches 
propagate plays the role of time and the perpendicular 
coordinates the role of space. Our 2D unloading sandbox 
is equivalent to a ID growing interface. Such an inter- 
pretation makes sense only when the stability condition 
and the avalanche dynamic rule is directional and local in 
space-time, such that causality is not violated in the in- 
terface growth interpretation. The stability condition (|l|) 
and toppling rule (|2|) of our model are row-by-row in na- 
ture and therefore indeed Markovian in this sense. 

Every stable configuration of the sandpile represents 
a possible interface growth life line (space-time-evolution 
interface world sheet). The conventional procedure for 
determining the scaling properties of growing interfaces 
is to average over a large set of completely independent 
MC runs. This would mean, in sandbox language, an en- 
semble average over completely refreshed surfaces, each 
totally uncorrelated from the previous one (except typi- 
cally for the initial condition in row y = 0). The toppling 
rule (||) is applied to all sites in every row, and repeated 
row-by-row, instead of only the unstable sites created by 



4 



time=t 

| h - - - - min[h(x- 1 ),h(x+ 1 )] 




x 



FIG. 3: The interface growth dynamics described by Eq. (y) 
with upper panel showing movement of steps (from the drawn 
to dashed line) and lower panel random depositions (shaded 
area) to the interface 

toppling only the highest site in the initial row. 

In avalanche dynamics, however, two subsequent grow- 
ing interface life lines in this ensemble differ only inside 
the avalanche area. From the interface growth perspec- 
tive this represents a rather peculiar and dangerous cor- 
related type MC run averaging procedure. The MC runs 
of KPZ space-time configuration are strongly correlated, 
and this raises the specter of a change in the interface 
roughness scaling properties. The numerical evidence, 
presented below is sufficiently ambiguous that this issue 
will preoccupy us in the second half of this paper. 

IV. KPZ GROWTH 

In this section we demonstrate that the interface 
growth model conjugate to the unloading 2D sandbox 
belongs to the 1+1D KPZ universality class. The time 
evolution of the interface is governed by the the toppling 
rule of the sand model with y in Eq. (Q) representing time 
t, 

h(x,t + 1) = nan[h(x + l,t), h(x - l,t)] + t](x,t). (3) 

In the conventional global type interface evolution (i.e., 
totally refreshing non-avalanche-type uncorrelated MC 
runs) every site in row t + 1 is updated according to this 
rule. 

Figure || illustrates the interface dynamics for one time 
step, t — » t + 1 . Conceptually the time step can be split 
into two parts; the deterministic min[ ] operator part and 
the stochastic random deposition r\ part. 

Note that because of the diagonal orientation of the 
square lattice (see Fig. ||), the lattice sites are not "sta- 
tionary in time" . The conceptually easiest interpretation 



to resolve this flip-flopping is to first double the number 
of lattice sites and then to require them to be paired al- 
ternately with their right or left neighbors at even and 
odd times; at even times sites 2n and 2n + 1 are fused to 
be at equal heights and at odd times the 2n — 1 and 2n 
sites. 

The upper panel shows the deterministic first half of 
the update (from the drawn to the dash line). The 
partners switch and the min[ ] operation equalizes their 
heights by choosing the lowest of the two. so this step 
always removes material. 

This can be interpreted also in terms of a movement of 
the steps in the interface. All up-steps move to the right 
and all down-steps to the left; while up and down steps 
merge when they meet at one site. 

The lower panel illustrates the second half of the up- 
date. The height of each fused pair increases by a random 
amount < r\ < s c . 

Deposition-type interface dynamics like this typically 
belongs to the KPZ universality class fl^| . Indeed, 
Eq. (^) can be rewritten as 

h(x,t + l) = -[h(x+l,t) + h(x-l,t)\ 

- - \h(x + l,t) + h{x - 1, t)\ + r)(x, t),(4) 

and from this easily be identified to be a discrete form of 
the KPZ Langevin equation, 

f t ^h-\{Vhf + v . (5) 

The crucial point is that the coefficient of the nonlinear 
term A is clearly present. There is no hidden special 
symmetry of some kind that makes it vanish by accident. 
At A = 0, the KPZ equation would reduce to EW growth. 

To confirm the KPZ nature and make sure that the 
A is large enough that corrections to scaling from the 
EW point (A = 0) are not obscuring the KPZ scaling, 
we perform MC simulations on the interface dynamics 
as illustrated in Fig. |3[ The MC runs are completely 
independent. 

We measure the time evolution of the interface width 
W defined as 

W 2 (L x ,t) = ((h-h) 2 ) (6) 

with over bars (angle brackets) indicating average over x 
(ensemble). Starting from, e.g., a flat initial condition it 
should scale as 

W ~ iP (7) 

at intermediate times <C t <C L x , and saturate at 

W ~ L% (8) 

for t 3> L x ; with L x the length of the ID interface. The 
exponents for the KPZ universality class in 1+1 D are 
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FIG. 4: MC results for the global interface width: left, finite- 
size (L x ) estimates for the saturated surface width expo- 
nent a; right, finite-time estimates for the transient interface 
width exponent (3 from a flat initial configuration. The solid 
(dashed) curves are for continuous (discrete) height model. 



known exactly with a — 1/2, [3 = 1/3, and z = a/f3 ~ 
3/2. 

The numerical results are shown in Fig. ^. The val- 
ues of ct{L x ) are obtained from the saturated interface 
widths by imposing the scaling form (||) at adjacent val- 
ues of the system size L x . Similarly, the values of (3{t) 
are obtained from the transient interface widths by im- 
posing the scaling form ((?]) at nearby times t. We like 
to remind the reader that simple log-log plots of W ver- 
sus L x and t look typically impressively straight but are 
notoriously inaccurate. The construction of effective ex- 
ponents, in the above manner might at first glance look 
less impressive (the data appears noisier) , but this brings 
the analysis to a higher level where the leading correc- 
tions to finite-size and finite-time scaling become visible. 

The approach to L x — > oo in Fig. ^ is consistent with 
the leading correction to scaling exponent y lr = —1/2 
expected from the EW term V 2 /i in Eq. (||). The correc- 
tions to FSS are stronger when the height variables are 
discrete than when they are continuous. This is consis- 
tent with the smaller growth rate in the discrete height 
interface, and the fact the growth rate is typically pro- 
portional to the nonlinear term A. On average, more 
material is removed during the first deterministic part of 
the update process when the surface heights are discrete. 



V. SCALING PROPERTIES OF 2D 
AVALANCHES 

In this section we derive the exact relations between 
the scaling properties of the avalanches and 1+1D KPZ 
interface growth. However, in the latter the world sheets 
are sampled in the correlated manner as outlined in 
Sec. flU 



The characteristic feature of SOC is the lack of typi- 
cal avalanche length, width, depth, or mass scales. The 
probability distributions follow power laws. For example, 
the distribution of avalanche widths scales as 



with scaling exponent r w . Similarly, the avalanche 
length, depth, and mass distributions scale as power laws 
with exponents 77, t$ and r m . We can summarize this in 
a metadistribution function P(Z, w, 5); the probability to 
find an avalanche of a specific width w, length I, and 
depth S, obeys the scaling relation 



P(l, w, 5) = b- a P{b- z l, b^w, b- a 6) 



(10) 



with b an arbitrary scale parameter. The exponents a, z, 
and a are expected to be robust with respect to details 
of the dynamic rule, and thus are characteristic of the 
universality class to which this avalanche dynamics be- 
longs. Single parameter distributions, such as P w , follow 
by integrating out the other variables. This implies the 
following expressions for the r exponents, 
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or inverted, 
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(11) 



(12) 



Let's presume that the avalanches are compact, i.e., 
that the inside and the boundaries of an avalanche are 
well defined and distinguishable (unlike in certain frac- 
tal structures), and that the sizes of the holes (unaf- 
fected regions) inside the avalanche do not scale with the 
avalanche size. This can be checked visually from typi- 
cal simulation configurations, and both assumptions are 
indeed satisfied in our dynamics at least qualitatively. 
In that case, the mass of the avalanche must scale as 
to ~ IwS, such that the critical exponent of the distribu- 
tion of avalanche masses P m ~ mT Tm obeys the identity 

(13) 



1 + z + a 



There is one more relation between these critical ex- 
ponents (leaving only two independent ones). The 
avalanche is initiated by lowering the bar at the driv- 
ing edge of the box. In the stationary state the average 
surface profile is invariant, and therefore it shifts down 
at the same rate as the lowering bar. Thus we know how 
much mass drops out of the box on average. 

To be more precise, during each avalanche event, the 
height of only one single boundary site at y = is low- 
ered by, on average, an amount s c /2. For a sandbox of 
width L x the boundary row is lowered by s c /2 after L x 
avalanches. In the stationary state, the entire surface 
matches this lowering speed, such that the amount of re- 
moved sand is on average equal to L x L y s c /2. Therefore, 
the average mass of each avalanche must be equal to 



(m) = -s c L y . 



(14) 



The scaling properties of the mass distribution function 
tie into this because 



p 



w 



(9) 



(to) = / to' P m (m')dm' , 



(15) 
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which can be evaluated using the metadistribution func- 
tion as 

rL y />oo />oo 

(m) ~ / dl I dw I dSlwSP(l,w,S) 
Jo Jo Jo 

/>oc />oo />oo 

+ m Ly \ dl dw dSP(l,w,5). (16) 

V J Ly JO JO 

This equation incorporates finite-size effects. The box 
is presumed to be wide and deep enough, such that the 
the length L y of the box (in the direction perpendicu- 
lar to the driving edge) is the only limiting finite-size 
factor. The first term in the above equation accounts 
for all avalanches that fit inside the box and the sec- 
ond term for the ones that reach the L y edge, and thus 
are prematurely terminated. The first integral scales as 
Ly a+2+2z+2a )/ z f or large L y . The second term scales 
with the same power, because the second integral scales 
as L y a + 1 + z + a )/ z wn ile the mass factor in front of it 
scales asm~ IwS ~ L y 1+z+a ^ z . The result 

(m) ^ £(-«H-2+2*+2a)/» j ( 17) 

when compared to Eq. (|l4|), yields the exponent identity 

a = 2 + z + 2a. (18) 

The validity of these exponent identities goes well be- 
yond our KPZ type unloading sandbox. For example, 
the EW type directed avalanche models by Paczuski and 
Bassler M and Kloster et al. |l(| obey our Eq. ([ll]) when 
we substitute for z and a the EW values (z = 2, a = 1/2). 
The scaling exponents of the original Dhar-Ramaswamy 
model can be described by the same equations with z = 2, 
a = as well. 



VI. NUMERICAL RESULTS FOR 2D SANDBOX 
AVALANCHES 

The discussion of the previous section leaves us with 
two independent avalanche critical exponents, a and z. 
The notation anticipates their identification with the 
scaling properties of a rough interface in interface growth. 
There, a is the scaling exponent of the interface width 
and z the dynamic critical exponent. Indeed, the inter- 
face width relates to the depth of the avalanche, and time 
to the the length of the avalanche. We expect therefore 
that a and z take same values as in 1+1D KPZ growth, 
a + z = 2 and a = 1/2. 

We perform MC simulations on the sandbox avalanche 
model and measure the avalanche metadistribution func- 
tion P(l, w, S\L y ), see Eq. (|l0[). The sandbox is always 
taken wide and deep enough such that the box length L y 
acts as the only FSS type limiting factor. We average 
over 2 31 avalanches. The reduced distributions, such as 
Pi ~ l~ n , follow from the metadistribution from, e.g., 
summation over w and 5. 




FIG. 5: FSS plots for the r exponents of 2D sandbox model. 
The solid (dashed) lines are for continuous (discrete) height 
model. 

Figure ^ shows FSS approximates for the r exponents. 
They are constructed as follows. Power-law-decaying 
objects such as P/ ~ l~ Tl are almost always subject 
to crossover-scaling- type effects, i.e., subdominant addi- 
tional power-law terms. In the language of renormaliza- 
tion theory they originate from so-called irrelevant scal- 
ing fields and also from nonlinear scaling field effects. 
This is well documented in equilibrium critical phenom- 
ena, but most recent nonequilibrium scaling studies ig- 
nore this systematic effect, e.g., by simply making a log- 
log plot of Pi as function of I and drawing a least-square- 
fitting-type straight line through the data. Such results 
show very little statistical noise, but can give rise to sig- 
nificant systematic errors. An example of the importance 
of corrections to scaling, was the large spread in reported 
values of the stationary state roughness exponent a be- 
tween various 2D KPZ-type-growth lattice models, which 
was resolved using a similar FSS analysis as presented 
here Q. 

In the limit of large I the subdominant additional 
power-law terms fade away. So, more weight must be 
put on the large / part of the data than on the short I 
section. However, it is a balancing act, because at large 
I the results become noisier, since few avalanches reach 
that far. 

The total number of avalanches that reach beyond y 
scales as 

Q l (y)= / m)dl ~ - y-T>+\ (19) 
Jy T i 

if the fraction of avalanches of length y scales as Pi ~ 
A y~ n (these are only the leading terms). We construct 
a y dependent approximate for the exponent tj from the 
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FIG. 6: Effective scaling exponents derived from station- 
ary avalanche distributions of sandbox systems. The solid 
(dashed) lines are for continuous (discrete) height model. 



of the deviations makes the latter more likely (except 
when this happens to be a continuously varying expo- 
nents scenario). 

We will blame the correlated MC averaging feature for 
this, but it should be noted that avalanche distributions 
are intrinsically more sensitive to FSS effects than global 
interface features. Many avalanches in the ensemble are 
small compared to the global box size, and therefore sam- 
ple and average the KPZ scaling properties over much 
smaller lengths and shorter time scales than in a conven- 
tional global interface roughness analysis at a comparable 
space-time box size. 

One option is to push the run button on the computer 
and out perform all corrections to FSS. Unfortunately it 
would require extremely long MC times, to create large 
numbers of such large avalanches. It is doubtful we would 
be able to get far enough in a reasonable time span. 
Moreover this approach is intellectually unappealing. We 
prefer to search for the origin of the deviations in the ex- 
ponents. 



VII. AVALANCHE CORRELATED MC RUNS 



ratio of these two quantities, as 



n(y) 



iPijy) 
Qi(y) 



(20) 



The results are shown in Fig. |. (We do the same for the 
other distributions.) Plots such as this are intrinsically 
noisier than conventional simple log-log type of plots of 
the distributions, but they contain much more informa- 
tion. The variation with y reflects the leading correc- 
tions to scaling. The statistical noise at large y could be 
suppressed by running the MC simulation longer. The 
simulation time is the only limiting factor. We used 2 31 
avalanches and in that case, L y = 512 is the optimal box 
size. 

In Fig. [| we replot the same data in terms of a, z, 
and er, following Eq. ([l2]) and using the same type of 
FSS analysis. From the trend of the curves, we conclude 
that a = 0.46 ± 0.01, z = 1.52 ± 0.02, a = 4.43 ± 0.05, 
and T m = 1.48±0.01. This means that the exponent rela- 
tions ([[3]) and ([l8]) are satisfied well within the statistical 
noise limitations, i.e., within a few percent. 

Surprisingly, the actual values for z and a, although 
close, differ significantly from the exactly known 1+1D 
KPZ values, a = 1/2 and z — 3/2. They deviate more 
than warranted from statistical noise alone, and do not 
converge smoothly if the KPZ values are correct. The 
approximates for a actually undershoot the KPZ value 
a = 1/2, and those for z overshoot z = 3/2. This sys- 
tematic effect needs to be explained. It could be that 
the exponents differ in a fundamental manner from the 
conventional KPZ values, or that we are looking at un- 
usually large and slow corrections to FSS. The smallness 



The basic premise of our exponent identities is that 
avalanches are like any other fluctuation on a 1+ ID KPZ 
type world sheet. Initially flat KPZ interfaces (the sand 
surface next to the driving edge) roughen in time (mov- 
ing away from the driving edge) in such a manner that 
at (KPZ) time y the stationary state roughness is estab- 
lished within a length scale l x ~ y 1 /". This defines a 
so-called spreading cone. The avalanches are expected 
to follow the same pattern. However, the avalanche cone 
seems to spread slightly faster, since the above avalanche 
value for z slightly exceeds the conventional KPZ value, 
and inside the avalanche the surface seems to be slightly 
less rough, since the avalanche value for a is slightly 
smaller. 

In this and the following section we will establish that 
this is caused by correlations with previous avalanches. 
The new avalanche does not run its course on a pris- 
tine fresh KPZ interface world sheet but on an aged one 
scarred by previous avalanches. 

There are two obvious tests to address the effects of 
these scars. The first one is to determine the avalanche 
distributions for only the first avalanche on a fresh KPZ 
world sheet (the initial condition), i.e., to refresh the en- 
tire surface completely after each avalanche. The results 
are shown in Fig. (?[ The first-avalanches likely follow nor- 
mal KPZ exponents: z converges now smoothly towards 
z = 3/2; while the FSS approximates for a, although 
still too small, start to turn towards a — 1/2, and do 
not cross that value anymore. It should be noted that 
the FSS corrections are expected to be larger, and that 
the data is noisier than in Fig. ||, because although we 
ran the same number of avalanches (2 31 ), the fraction of 
large avalanches is smaller, leading to smaller and noisier 
amplitudes in the power-law tails of the distributions. 
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FIG. 7: Effective scaling exponents derived from the distri- 
butions of first avalanches on fresh sandbox surface for the 
continuous height model 



The second test of the role of the scars is to measure the 
global interface roughness for avalanche type correlated 
MC runs instead of completely refreshing MC runs. The 
upper panel of Fig. || shows the global interface width W 2 
as function of time for several L x s. The drawn lines cor- 
respond to avalanche correlated MC runs and the dashed 
line to conventional uncorrelated MC averaging. The 
drawn lines have bumps, i.e., the avalanche correlated 
runs lead to rougher interfaces at intermediate times. 

This enhanced interface roughness is caused by the 
scars left by earlier avalanches. The scars vanish at very 
large y because avalanches reaching that far span the en- 
tire system in the ^-direction. Figure |^ shows a typical 
configuration of scars. The lines are the traces of previ- 
ous avalanches, i.e., their edges. Latter avalanches wipe 
them out partially. 

For finite system sizes, the stationary state interface 
width follows from the plateaus at large times. There the 
avalanche correlated and uncorrelated MC curves coin- 
cide. This is to be expected, because the large avalanches 
that span the entire system (in the x direction at large y) 
occur at regular MC time intervals, such that the large y 
part of the surface (i.e., the stationary state of the growth 
process) is completely refreshed periodically and there- 
fore sampled effectively like in uncorrelated MC runs. As 
a result, the roughness exponent a, defined by Eq. (||), 
is the same for the both cases. 

Most avalanches do not extend into that large y part 
of the surface. They terminate in the scarred part of the 
surface. Therefore, we define an alternative roughness 
exponent a*, associated with the scaling of the bumps, 
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FIG. 8: Upper panel: square interface width for stationary 
sandbox surface (solid lines) comparing with fresh surface 
(dashed lines); Lower panel: the difference between the two. 
From bottom up, the corresponding system sizes, L x , in the 
transverse (x) direction are 8, 16, 32, 64, 128 and oo. 




FIG. 9: A typical configuration of the scars on the sandbox 
created by the avalanches. The driving edge is located at the 
bottom of the graph while avalanches propagate upward in 
the y (or t) direction. The system sizes are L x = 256 and 
L y = 512. 

in terms of the maximized width 

W* = max W(L X , y) ~ L x ' (21) 
y 

more relevant for the avalanche scaling properties. Note 
that for uncorrelated MC runs, a* = a, since the inter- 
face width increases monotonically in time. 

The conventional method for measuring the exponent 
(3, involves the slope at times y < L z x , and thus is sensi- 
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0.45 




FIG. 10: Finite-size approximates of the scaling exponents 
for stationary surface of sandbox (or correlated MC runs for 
the interface model) with a* defined by Eq. (|2l]) and /3 by 
Eq. (|). The solid (dashed) curves are for the continuous 
(discrete) height model. 



tive to the bumps in W as well. The results are shown 
in Fig. [l0| Compared to those in Fig. ||, they clearly 
converge less smoothly, with larger corrections to scaling 
and we should wonder if they converge to the conven- 
tional exact KPZ values, a = 1/2 and (3 — 1/3, at all. 

In the lower panel of Fig. @ we plot AW 2 as func- 
tion of time, the difference between the squared widths 
of avalanched correlated MC runs (the drawn lines in the 
upper panel) and completely uncorrelated MC runs (the 
dashed lines in the upper panel). For infinite system size, 
AW 2 scales as AW 2 ~ y s with an exponent that numer- 
ically is very close to s ~ 1/3. Since the width itself 
scales as W 2 ~ y 2 ^ 3 , it follows that the bumps in the 
width curves are a transient FSS effect. 

This settles our basic issue at the numerical level; the 
avalanche correlated nature of the MC runs does not 
change the interface scaling exponents, but only gives 
rise to slow corrections to FSS. In the next two sections 
we will identify these corrections to scaling with the scars 
on the surface left behind by previous avalanches. 

We start this analysis here by casting the deviations 
into the framework of corrections to scaling from a so- 
called irrelevant operator in the sense of renormaliza- 
tion theory. Let O sc (x) be that irrelevant operator and 
u be its scaling field. This mounts to presuming that 
the avalanche correlation between MC runs can be rep- 
resented effectively by adding to the KPZ Langevin equa- 
tion (|5|) , a term uO sc (x) . We will have to determine below 
how O sc (x) is related to the density of scars on the inter- 
face space-time world sheet left by previous avalanches. 
According to scaling theory, the presence of such a term 
to the Langevin equation leads to corrections to scaling 
in the interface width as 

W 2 (L x ,y, u) = b 2a W{b~ l L x , b~ z y, b y -u), (22) 

co, to 

W 2 (y,u) = y 2a ^S(yy^ z u) 1 (23) 

and by expanding the scaling function S, while assuming 
that y sc < 0, such that u = is a stable fixed point, and 



i.e., in the infinite-size limit, L x 



the argument y Va °/ z u is a small parameter, to 

W 2 (y, u) = y 2a/z [5(0) + y v ^ /z uS'{0) + •••]. (24) 

The critical exponent y sc of this irrelevant scaling field 
must take the value y sc — —a to account for the 5W 2 ~ 
yi/3 corrections in the interface width we found above. 
Moreover the operator must scale as 



O sc {x) 



(25) 



with critical dimension x sc — z, since the KPZ equa- 
tion (^), implies that the terms uO sc {x) and dh/dt must 
scale alike. In the following two sections we will trace 
down the geometric identity of this mysterious operator 
O sc , starting with the ID version of the model. 



VIII. SURFACE ROUNDING IN THE ID 
UNLOADING SANDBOX 

The ID version of the unloading sandbox shows 
the same type of differences between uncorrelated and 
avalanche-type correlated MC runs as the 2D version. 
We determined numerically the difference between the 
interface width for avalanche-correlated and uncorrelated 
MC runs, and found that it diverges as a power law 
SW 2 ~ y 1 / 2 , with an exponent which is again (like in 
2D) half the size of that for W 2 ~ y itself. According the 
corrections to scaling formalism (24), the scaling dimen- 
sion of O sc must therefore be equal to x sc = z, just as in 
2D. 

The underlying interface dynamics becomes a zero di- 
mensional growth model, i.e., a simple random walk in 
the h direction with a nonzero drift velocity to account 
for the net tilt of the surface. The exponents of the vari- 
ous avalanche distribution functions must obey the same 
type of relations as in Sec. 



a 



n 



a + z 



and 



2a 



(26) 



(27) 



Without loss of generality we can set a = 1 (measure 
all lengths in terms of 6). These identities are satisfied 
exactly, and the exponents are the same for uncorrelated 
and avalanche correlated runs. From the interface dy- 
namics perspective, a single directed random walker, the 
diffusion equation character of the dynamics implies that 
z = 2a = 2. The values of all the other exponents fol- 
low from this, and are consistent with their values from 
the avalanche perspective. There, we are dealing with 
the statistics of merging random walkers. The number 
of walkers at a given "time" y is equal to the number of 
avalanches of a length I equal or larger than y in the en- 
semble of MC runs. The density of the walkers decays as 
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p(y) ~ V such that the distribution of avalanche 

lengths obeys the form 



W) 



d_ 

dy 



p(y) 



r 3 / 2 , 



(28) 



y=i 



and therefore that 77 = 3/2. The depth of the avalanche 
follows from the maximum separation between two sub- 
sequent walkers, and scales as S ~ I 1 / 2 , i.e., a/z = 1/2. 
The mass scales as m ~ £5 ~ Z 3 / 2 , i.e., (a + z)/z = 3/2 
and r m = 4/3. 

This can be compared directly with the exponents of 
other ID sandpile models, e.g., with results by Paczuski 
and Boettcher Q on the so-called Oslo sandpile model 
where r = T m ~ 1.55 and D = (a + z)/z w 2.23. 

Let's turn our attention now to the central issue, 
the difference between uncorrelated versus avalanche- 
correlated MC runs. Adding a term like uO sc to the 
diffusion equation of motion creates a correction to the 
drift velocity of the random walk. This suggests we can 
identify the geometric meaning of O sc directly by study- 
ing the deviations of the slope near the driving edge of 
the surface from its asymptotic value. 

The average surface slope does not show any devia- 
tions (near the driving edge) from s c /2 when we run the 
dynamics as a conventional random walk, which amounts 
to "completely refreshing" the surface after each MC run 
(uncorrelated MC runs). The avalanche-correlated runs 
do show a rounding of the surface near the driving edge 



s(y)~Ay K + 



1 

2 



(29) 



The numerical results for the exponent yield k = 0.98 ± 
0.03, in accordance with x sc = z from the interface width 
since k — x sc /z and z — 2 for random walks. 

This rounding originates from the distribution of ter- 
mination points of the avalanches. A new random walk 
starts below the previous one and propagates until it 
meets the previous trajectory and terminates. The 
avalanche is the space between the trajectory of that 
new random walk and the already existent surface. The 
amount of rounding of the slope near the driving edge is 
proportional to the distribution p(y) of merging points 
on the surface. Those are the scars from previous 
avalanches. Each random walker by itself does not con- 
tribute to the rounding, i.e., on average its walk has 
constant slope s c /2 independent of time. However the 
merging process truncates each walk and does so in an 
upwards biased fashion. Each merging event causes the 
surface to drift upwards by a certain amount (s c /2, in av- 
erage, for the discrete h version). Therefore the rounding 
of the surface is proportional to p{y) . 

The entire process and the set of subsequent stable 
sand surfaces (Fig. |ll|) is therefore equivalent to a system 
of merging random walkers obeying the rule A + A — > A. 
That type of dynamics has received extensive attention 
recently and its various scaling properties are known ex- 
actly (nJ. There is little doubt that our ID unloading 




FIG. 11: Traces of stable sand surface over 256 avalanches for 
ID sandbox model with L y = 256. The system is driven from 
the left at y — 0. 



sandbox is exactly soluble, using absorbing wall type ran- 
dom walk mathematics |17j. However, we will refrain 
from pursuing this path in this paper. 

The critical dimension of O sc ~ p(y) can be estimated 
(for intuition building purposes) as follows. After adding 
a term uO sc to the KPZ equation we should also write 
down an equation of motion for O sc itself, to close the 
equations. The latter is not trivial, because the scars on 
the surface build up slowly in time, such that that the 
equation of motion for O sc is highly nonlocal. On the 
other hand, the linear nature of the diffusion equation 
allows one to be somewhat frivolous with the order in 
which averages are taken, (without losing the essential 
physics, nor even the correct critical exponents). 

Let Pt(y) be the endpoint distribution after t 
avalanches (MC time steps). During the last MC time 
step, one avalanche runs through the system. It refreshes 
the entire surface before its termination point y = l~ t . 
such that pi at site y does not change if the avalanche 
terminates before y; pi(y) = 1 if it terminates at y; and 
Ptiv) = if it extends beyond y: 



dPiilj) 
dt 



Pi(y)-P~ t (y) W)dl 



(30) 



with Pi{l) the probability that the avalanche terminates 
at distance / from the driving edge. The stationary state 
endpoint profile therefore takes the form 



p{y) 

and Pi{l) ~ l- Tl yields 

p{y) 



Pi{y) 



1 



-y 



(31) 



(32) 



n-r 

In other words, the surface curvature scales as As ~ 
y-x sc /z Xsc _ Z: in agreement with the above re- 
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FIG. 12: Scaling exponent for boundary correction to the 
local slope of fresh 2D sandbox surface (or, in the interface 
language, transient growth rate from a flat interface), S{(y) — 
Sf (oo) ~ y~ Kt , and its correction due to the iterated avalanche 
process, As = s(y) — Sf(y) ~ y~ K '. 



suits. Interestingly, this result is independent of the 
actual value of the scaling exponent 77, provided that 
t; > 1, which has to be true for Pi to be normalizable. 

In conclusion, in ID we identified the crossover scaling 
operator with the density of avalanche endpoints. These 
represent indeed the scars on the surface, the memory of 
previous avalanches. 



IX. AVALANCHE ROUNDING NEAR THE 
DRIVING EDGE IN 2D 

As in the ID model, the surface slope is modified by 
the iterated avalanche process. However, unlike in ID, 
the average slope near the edge is not constant already 
in conventional interface dynamics (where the entire sur- 
face is being refreshed during each MC run) . The surface 
slope is related to the growth rate of the underlying inter- 
face, and the rounding of the slope near the driving edge 
represents the transient growth rate of the KPZ interface 
from the initial configuration, e.g., a flat one: 



Sf{y) ^v a + cy 



(33) 



with y playing the role of time and the subscript, f, de- 
noting that the entire surface is refreshed. By direct nu- 
merical simulation of uncorrelated interface dynamics we 
find K[ « 0.7 (the left panel of Fig. |l2|). This is consistent 
with conventional KPZ scaling and power counting 



h/y ~ y a ' z 



I) 



-2/3 



(34) 



suggesting m — 2/3. 

We evaluate the surface slope profile s(y) in avalanche 
correlated dynamics MC runs, in terms of the difference 
with respect to the uncorrelated case, 



As(y) = s(y) - s { (y) ~ y H 



(35) 



The FSS analysis for the exponent n (the right panel of 
Fig. ^ yields K = 1.05 ± 0.07. This is in agreement 
with .t sc = z and k = x sc /z implied by the corrections to 
scaling formalism (pi|). 



Old 
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FIG. 13: Two possible cases at a boundary of an avalanche 
cluster (the shaded area) : a. avalanche expands; b. avalanche 
shrinks. The local slopes along the arrow marks is reduced in 
a. while increased in b. 



Inside the bulk of an avalanche the interface is fully re- 
freshed, and scales as in uncorrelated KPZ dynamics. At 
the avalanche boundaries, the slope of the surface is bi- 
ased upwards, because of the merging with previous MC 
runs (which are on average shifted upwards by an amount 
s c /2L x each time an avalanche is triggered). This means 
that the As is proportional to the density of scars in the 
surface. In ID, the scars are point-like objects, the end- 
points of the avalanches; but in 2D the avalanche bound- 
aries are line objects. This nonscalar aspect makes that 
most line-segment contributions, when integrated along 
the boundaries of an avalanche, cancel out against each 
other. 

To be more precise, s(y) represents only the compo- 
nent of the slope in the y-direction, and the magnitude 
of those jumps depends on the local angle 9 the bound- 
ary makes with the y-axis. This is an odd function, 
A(0) = -A(-0), as illustrated in Fig. [H| The slope 
change is negative when the avalanche opens up and pos- 
itive when it narrows down. The latter also implies that 
A(0) has opposite sign for the left and right boundary of 
each avalanche. Notice that, while in the lattice model 
takes only two discrete values, it renormalizes to a con- 
tinuous variable at larger length scales. 

Let's estimate the change in surface slope due to these 
scars in the same spirit as we did successfully in ID. 
Consider one specific surface, and let sj(y) be the sur- 
face slope in a slice of the surface at distance y from the 
driving edge, averaged over all x, after t avalanches (MC 
time i). The last avalanche changes this as follows. Let 
wi(y') be the width of this avalanche, which terminates 
at y = If, in slice y' . The inside area of the avalanche 
is completely refreshed and therefore has the same av- 
erage slope Sf(y) as in ordinary KPZ dynamics (totally 
refreshed subsequent world sheets). This leads to the 
following equation of motion, 



dt 



[A(0 L ) - A(0 R )] + wiivftativ) - Si (y)} (36) 



The first term on the right hand side represents the 
creation of the two new avalanche edges, and the sec- 
ond term represents the refreshed surface inside the new 
avalanche. Note that ds^(y)/di = when this latest 
avalanche does not reach slice y, and that this is auto- 
matically taken care of because in that case 9^—9^ — 
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and A(0) = 0, while w^(y) = for y > l~ t . In the sta- 
tionary state, after averaging over all possible avalanches, 
Eq. (pq) leads to 



wdv) [* (v) - s-M] = a(0l) - A(e R ) (37) 

Next, we perform an heuristic coarse-graining 
renormalization-type transformation. At large length 
scales, the average angle 9 remains small, such that the 
right hand side can approximated as 



A(6 L ) - A(6 R ) ^ae L -6j 



dw~ t {y) 
dy 



(38) 



Finally, we presume that in the stationary state it is not 
too bad to treat the KPZ height fluctuations deep in- 
side the bulk of an avalanche and those near its edge as 
decoupled (at least in lowest order) such that 



d 



As(y) = s f (y) - s t (y) = a — log(w t ~(y)). 

dy 



(39) 



This yields As(y) ~ y^ 1 , exactly the power-law decay we 
are looking for, and consistent with all the above numer- 
ical results. 



The only requirement for the latter is that u>f (y) ~ y 



decays as a power law. Again, like in Eq. (J32j) for ID, the 
value the critical exponent £ does not matter. w^(y) is 
equal to the average avalanche width in slice y averaged 
over all avalanches. It is reasonable to expect, and we 
confirmed numerically, that this quantity scales with the 
same exponent as the average width of all avalanches 
longer than y, i.e., as 



w(l)P{l)dl ~ y 



l/z-n + l 



(40) 



which yields £ ~ 1/3. 

We are now ready to represent the crossover scaling op- 
erator O sc (x) in terms of the scars on the surface. Con- 
sider time slice y. O sc (x) = when no scar line runs 
through site x, and otherwise is proportional to the angle 
the scar line makes with respect to the y-axis. However 
the sign also flips depending on whether this represents a 
left or right boundary of the original avalanche. The lat- 
ter can be denoted by an arrow along the avalanche scar 
line. Alternatively, we can associate an age-field g(x,y) 
to the entire surface, representing the age of the surface 
segments (how many MC time steps ago site x was up- 
dated), 



|Vff| 



(41) 



with i y a unit vector in the y-direction. The denominator 
arises because the magnitude of the age jump across the 
scar line |Vg| docs not play a role. 



X. SUMMARY 



In this paper, we studied a directed avalanche model 
inspired by the unloading of a sandbox by means of a 
slowly lowering wall, and the wish to setup an avalanche 
dynamic rule belonging to the same universality class as 
KPZ type interface growth. The 2D sand surface repre- 
sents the world sheet of the 1+1D growing interface. 

The scaling exponents of the avalanche distributions 
are directly related to the dynamical and stationary state 
roughness exponents z and a of KPZ growth in 1+1D, 
Eq. (11). However, we encounter one crucial difference. 
From the avalanche perspective the conventional uncor- 
related MC runs correspond to completely refreshing the 
surface, i.e., an ensemble average over all possible initial 
conditions, without ever running an avalanche. From the 
KPZ perspective, the avalanche dynamics represents an 
unusual MC ensemble averaging procedure where subse- 
quent interface world sheets only differ inside the single 
avalanche. This avalanche-correlated-type averaging en- 
hances the interface roughness at time scales y < L*, due 
to the scars of previous avalanches. It required a care- 
ful study, combining numerical and analytical tools, pre- 
sented in the second half of this paper, to establish that 
these scars give rise only to larger than usual corrections 
to scaling and not to fundamentally different values of 
the global roughness scaling exponents z and a. 

The effect of the scars can be represented by introduc- 
ing an additional age field g[x,y) to the height variables 
h(x,y), that keeps track of how many MC runs ago site 
(x, y) participated in an avalanche. This age-field cou- 
ples into the KPZ equation (|J) as an additional term of 
the form uO sc . The operator O sc is proportional to the 
angle a scar makes with respect to the time-axis, and 
can be expressed in terms of the age field as shown in 
Eq. ( fll| ) . We establish that the coupling of this age field 
to the KPZ equation is irrelevant in the sense of renor- 
malization theory, both numerically and by writing down 
approximate equations of motion for uO sc . The scaling 
field u renormalizes with exponent y sc = —a and O sc 
scales with critical dimension x sc = — z. 

We believe that the results of our work presented here 
can be generalized to most "Markovian" avalanche dy- 
namic systems with local row- by-row type toppling rules, 
and that this is a promising route to improve our under- 
standing of the scaling properties of avalanche dynamics 
in general. 

This research is supported by the National Science 
Foundation under grant DMR-9985806. 
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